%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Explanation of what the file does
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This function solves for new revenues (hat or dots) in the RUV economy
% with forward looking migration using a linear system of equations
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Inputs needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Trade shares in levels (shareslev), size I*I*S
% 2. Final good demand (Vmat), size I
% 3. Final good consumption shares (Amatrix), size I*S
% 4. IO matrix (inoutmat), size I*S*S
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Outputs produced
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Revenue matrix (in levels, denoted Revmat), size I*N

function Revmat=NEBA_RevenueSolver(shareslev,Vmat,Amatrix,inoutmat)

% Define number of regions and sectors

I         = size(Amatrix,1);
S         = size(Amatrix,2);

% Revenues are linear in themselves, so I don't need a solver. However, the 
% matrices involved are complicated since some have three dimensions. 
% So I "vectorize" the (IxS) matrix of revenues into a (IxS)x1 vector of 
% revenues. Then, I can solve a linear system to obtain this vector. This 
% system will require to use the IxIxS shares matrix, so I make a huge 
% (IxS)x(IxS) matrix of shares that has in its S diagonal blocks the IxI 
% matrices for each of the sectors, that is the "shares big" matrix.

intermau  = shareslev;
qcell     = num2cell(intermau,[1,2]);
sharesbig = blkdiag(qcell{:});

% With this sharesbig matrix I can obtain the vector A and matrix B in the 
% system Revvec = A + B * Revvec, where Revvec is the revenue vector, A is 
% a (IxS)x1 vector, called "constvec" below, and B is a big (IxS)x(IxS) 
% matrix called "coeffmat" below. I obtain Revvec from the system as 
% (I-B)^(-1)*A where I is a (IxS)x(IxS) identity matrix. I then turn the 
% (IxS)x1 revenue vector into a IxS revenue matrix

constvec  = sharesbig * (reshape(Amatrix,I*S,1).*repmat(Vmat,S,1));
lambdamat = reshape(permute(intermau,[2,1,3]),I,I*S)';
phismat   = reshape(permute(inoutmat,[2,1,3]),S,S*I);
coeffmat  = eye(I*S) - repmat(lambdamat,1,S).*kron(phismat,ones(I,1));
Revvec    = coeffmat\constvec;
Revmat    = reshape(Revvec,I,S);

end
